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O ! Abstract 
00 

We study a "burst" event, i. e. the evolution of an initial condition 
having support only in a finite interval of /c-space, in the continuum 
shell model due to Parisi p]. We show that the continuum equation 
without forcing or dissipation can be explicitly written in characteris- 
tic form and that the right and left moving parts can be solved exactly. 
| When this is supplemented by the appropriate shock condition it is 

' possible to find the asymptotic form of the burst. 
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1 Introduction 



In the study of fully developed turbulence the so-called shell models have 
become an important tool. In these models the essential ingredients are 
the conservation laws and the assumption of locality on a discrete set of 
exponentially growing momentum shells. The models are constructed such 
that Kolmogorov's 1941 scaling law is built in as a fixed point in the absence 
of forcing and viscosity. 

It has recently been found (see e.g. [Q] and references therein) that these 
models actually capture important features of the strongly intermittent be- 
haviour seen in experiments, which gives rise to corrections to the simple 
Kolmogorov dimensional scaling. In order to understand how such correc- 
tions come about, it is necessary to study the strong "burst" events which 
destroy the homogeneity of, say, the dissipation field and presumably lead to 
the appearance of multiscaling. 

In 1990 Parisi |1| suggested a continuum approximation of a shell model 
and showed that solitary wave excitations would exist on top of the Kol- 
mogorov spectrum (a short presentation of his approach can be found in 0). 
It has, however, recently been stressed that the shell models have another 
(trivial) fixed point, where the fluid is quiescent, and that this state seems 
to be approached before a large burst. In the present paper we shall thus 
study how a single pulse (a disturbance initially localized in wave number 
space) will propagate into a motionless fluid. The details of how such a burst 
propagates are quite interesting and unexpected and although it remains to 
be seen, whether the shocks which are generated represent the bursts in real 
turbulent fluids, we believe that they are worth discussing since the methods 
that we use might be generalized to more realistic models. We feel sure that 
Giovanni Paladin would have enjoyed our story - indeed he was the one who 
introduced us to the shell models and later contributed and inspired so much 
of the work in this field. 

2 The GOY shell model and its continuum 
limit 

The basic idea behind the shell model is that the Navier-Stokes equation, 
when considered in Fourier space, receives contributions to the time deriva- 
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tive of the velocity from the velocities in a triangle of fc-vectors. This con- 
dition is mimicked by taking fc-vectors as the one- dimensional discrete set 
("shells") k n = k r n , and introducing only a single complex field u n for each 
k with interactions only between nearest and next nearest neighbours. The 
field u n represents the characteristic velocity differences across the nth shell, 
and, in the limit of vanishing viscosity and forcing, energy conservation is 
imposed in the form 

\u n \ 2 = const (1) 

The most studied model is the so-called GOY model Q (after Gledzer, 
Okhitani and Yamada) 

,du n , o , . 5 1 — 5 , n , n , 

[—7T + vk n u n ) = -ik n (u n+1 u n+2 u n _iu n+1 — w„_iw n _ 2 ) + F n (2) 

at r r z 

where v is the viscosity and F n is the forcing. 5 is a free parameter, but 
becomes fixed at S = 1/2 when the conservation of helicity is also satisfied. 



2.1 The Parisi continuum limit 

It is an interesting question how to define the shell model in the limit where 
the fc-variable becomes continuous. A particular way of reaching this limit 
is to take r — > 1: By writing the distance between the shells as r = 1 + e 
with e C 1, we have k n ~ exp(ne), so with n ~ const. /e a continuous range 
of values is obtained for the variable k. 

To proceed we use a Taylor expansion of the type 

du 

omk 

ok 

= u(k,t) + ek — ' , with u(k, t) = u(k n , t) = u n (t), (3) 

OK 

and similarly for w n+2 ,-u n _i, and u n -i- To first order one obtains 

u\ + uk 2 u* = -ik{2 - 5) (u 2 + 3kuu k ) , (4) 

where u k = du/dk etc. The higher order corrections to @j and the conver- 
gence to the discrete model will be examined in H. By rescaling time with 



3 



2 — 5 we get the Parisi equation, which in its inviscid, unforced form is: 

u* t + 3ik 2 uuk = —iku 2 (5) 
In terms of the real and imaginary part (u = a + ib) it can be written as: 

a t — 3k 2 abk — 3k 2 ab k = 2kab (6) 

and 

b t -3k 2 aa k + 3k 2 bb k = k[a 2 -b 2 ). (7) 

For the case of a = 0, a number of exact solutions to (0) have been found in 
1- 

3 Solution of the inviscid Parisi equation 

We shall now study the inviscid, unforced Parisi equation (§) in the special 
case, where the initial condition is a single "burst", i.e. that the field u is 
only nonzero in some finite interval in /c-space. 

3.1 Direct simulation 

The Parisi-equation @ or is a hyperbolic equation ||, and hyperbolic 

equations are notoriously hard to solve numerically due to the appearance 
of shocks in the solution. These can be dealt with in a crude way by using 
upwind differencing of the flux terms j7| , but to do this, one need to figure out 
in which direction a shock is moving. For the real part (|6]) there is a minus 
sign on the flux term, signifying that a shock moves to the left (provided 
both real and imaginary parts are positive). The imaginary part (|7|) is not 
so simple, since the flux term is composed of both a right and a left moving 
part. Never the less upwind differencing works there too, by assuming that 
the dominant direction is to the right. 

An example of a numerical solution with semi-implicit update in time, is 
seen in figure |l|. The initial condition is the zero state, with a gaussian peak 
a(k,t = 0) = b(k,t = 0) = Cexp(-(A; - k ) 2 ), with k = 40 and C = 0.25. 
In order to neglect viscosity, the inertial subrange must be to the right of the 
initial perturbation. 
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Figure 1: A numerical solution of the Parisi equations, without forcing and 
viscosity. The upper figure shows the solution at t — 1CT 3 after having been 
started with a gaussian peak at k — 40. The inset on the uppermost figure 
shows the part moving to the left, and illustrate that the relation a = \^3b 
is fulfilled to a high accuracy. In the lower figure only the imaginary part is 
shown in a space-time plot. Here the splitting of the initial pulse into a left- 
and right-moving part is clearly visible. The initial phase is 7r/4. 
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A short while after the start, shocks are formed in both directions. As 
was seen from the equations, the real part is only able to move to the left, 
while the imaginary part moves both ways. Thus only the imaginary part is 
important in the inertial subrange. Setting a = in the original equations, 
we get just a single real equation for the inertial range: 

b t + 3k 2 bb k = -kb 2 (8) 

For the left-moving part of the pulse, the real and the imaginary parts are 
seen to become proportional. Inserting b = Ca into the Parisi equation gives 
us two possibilities: a = (which we already treated) and a = ±a/35. For 
a = a/36 the equation for the imaginary part becomes: 

b t - 6k 2 bb k = 2kb 2 . (9) 

The explanation for the splitting of the pulse in this way will be given in the 
next section by writing ([|) in characteristic form. 

Another problem occurring when solving hyperbolic differential equations 
is to find the correct position and velocity of the shocks. Strictly speaking, 
the equations are not well-defined after the occurrence of a shock, because 
the solutions become multivalued. To solve past the time when the shocks 
first appear one has two choices: 1) to derive a shock condition by imposing 
conserved quantities (note that energy conservation alone is not enough, since 
two conditions are needed) or 2) include a higher order terms to smooth out 
the shock. When solving differential equations with finite differences, there is 
always a small amount of "numerical diffusion", corresponding to a d 2 u/dk 2 
term in the equations. This is not necessarily the correct higher order term, 
and therefore the velocity of the shocks as seen in figure |T| is not necessarily 
correct. We shall show later how to obtain the correct shock conditions by 
imposing energy conservation, and thereby the correct asymptotics for the 
solution. 

3.2 Characteristic form for the Parisi equation 

Surprisingly, the "source" term (i.e. the right hand side) of equation @ can 
be removed by the simple substitution u = k~ x ^v and k = (2x)~ 3 ^ 2 . This 
turns (H) into a complex equation reminding strongly of the Burgers equation 

v * - ivv x = 0. (10) 
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In terms of real and imaginary parts of v (again denoted by a and b) we get 



a t + ba x + ab x = (11) 

and 

b t + aa x - bb x = (12) 

The calculations turns out to become simpler if (|T0|) is written in terms of 
modulus and phase, i.e. 

v = re** (13) 
whereby we find the equivalent equations 

ft + rr x sin 9 H — r 2 8 x cos 9 = (14) 
3 

and 

9 t + 3r x cos9 -r9 x sin9 = (15) 
where 9 = 30. This can be written as a matrix equation 

q t + Aq x . = (16) 

where q = (r, 8) and 

a _ ( rs i n ^ |?" 2 cos# \ , . 

I 3 cos 9 —r sin 9 I 

The characteristics are the eigenvalues of A, which are A = ±r and the 
Riemann invariants can be found through the left eigenvectors of A (or the 
right eigenvectors of the transpose of A [§]). These eigenvectors are of the 
form e = (a, (3), and satisfy 

/±l-sin#\ 
\ 3 cos 9 J 

The Riemann invariants satisfy 
and 

% = >> ™ 
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and the Parisi equation ([10]) can be reformulated as the two equations 

dx 

J± = const on the curve — = ±r(x, t) (21) 

Only the ratio {3/ a. is uniquely defined. One solution for ([0J)-(^(]) is obtained 
by taking a = 1/r and thus 

at = — ± = - 22 
ar r 

and 

_ ±l-sinfl 
^" 90 " 3cos£ ^ j 

which gives 

J ± = log(r(l±sinfl) 1/3 ) (24) 

Now, the Riemann invariants are only defined by their invariance along the 
characteristics, and thus we can equally well take them to be defined without 
the logarithm, i.e. 

J± = r(l±sinfl) 1/3 (25) 

which of course is constant on the same locus. 

To complete the solution we can express r and 9 through J±. We find 



(26) 



and 



si ^ = ls— (27) 



J%-Jl 
Jl + J* 

The formulation in terms of characteristics can now be used to understand 
the splitting of the pulse found in section 3.1. The important point is that, 
since r is nonnegative, one family of characteristics (J + ) can never move to 
the left while the other one (J_) can never move to the right. If the initial 
condition has support in a limited region of 1 the same is true 

of J±. They both vanish (in the initial condition) outside of this interval. 
For times t > we compute the field values by finding where, on the x-axis 
(i.e. t — 0) the characteristics going through the point (x, t) emanate. Now, 
if x > x + , the ./.-characteristic going through this point must emanate from 
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some xo > x + and thus J_(x, t) = J_(xo,0) = 0. This means that either 
r = (which makes the entire field vanish) or sin# = 1 which means that 
9 has the constant value tt/2 or = tt/6. This corresponds exactly to the 
result of last section: that a = \f?>b for the left moving pulse (which is moving 
right in the x- variable). 

A completely similar argument, valid for x < x_, shows that in this case 
sin# = — 1 or 9 = 3tt/2 and = 7r/2, which implies that a = 0, again in 
agreement with the result of last section for the right moving pulse. This 
shows that the regions outside are so-called simple wave regions [|]. 

When we said that 9 = tt/2 =^> = tt/6 it is not entirely correct. The 
value = 7i I 6 is only one possibility, the two others being = 5tt/6 and 
= 37r/2. Likewise, 9 = 3tt/2 can mean = tt/2, = 7tt/6 or = ll7r/6. 
If the initial phase is in the interval between two of these six values of 
corresponding, consecutively, to sin^ having the value +1 and —1, the values 
of selected will be precisely those two. The values = 7r/6 and = tt/2 
given above correspond to an initial phase around 7r/4. 

3.3 Asymptotic analysis of the burst 

The splitting of the pulse into a right-moving part with a = and a left- 
moving part with a = \/3& makes it possible to give a complete solution for a 
single burst event. We shall divide the analysis into the two simplified cases; 
forward and backward moving parts. In each regime the field is described 
by a single scalar equation. It is therefore possible to find the correct shock 
condition using a single conservation law: energy conservation. We have 
chosen the case where the initial phase is around 7r/4, but the other cases 
can be handled similarly. 

3.3.1 Forward propagation 

When the real part a = the equation for the imaginary part b (which we 
shall call u in this section) is 

u t + 3A; 2 uufc = —ku 2 (28) 

where we, for the sake of explicitness, use the original form instead of the 
"Burgers form" (|10|). Again, this is a hyperbolic equation, which can be 
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solved by the method of characteristics. The characteristic equations are 

k'(t) = 3k 2 u (29) 

u\t) = -ku 2 (30) 
and they can be solved e.g. if initial conditions are specified in the form 

u(k,t = 0) =uq(x) (31) 

where x = ko = k(t = 0) (note that this x is different from the one used in 
section 3.2). Indeed one can note that for U clS Sb function of k (which is what 
we want) (|29| - |30D leads to 



du u'(t) 1 u 



(32) 



dk k'(t) 3 k 
which can be integrated to: 

u = u {x)(k/xy 1/3 (33) 



This can now be used to solve (|29"| - |3~0D explicitly as: 

u{x, t) = Uq{x)(1 - 2f(x)t) 1/2 (34) 
k(x,t) =x{l-2f{x)ty 3 / 2 (35) 

where f(x) = xuq(x). 

Two types of singularities occur in these solutions: 

1. u(k) becomes multiple- valued when k'(x) = 0. This will be called a 
"turning point" and it implies the existence of a shock. 

2. Both k(x) and u(x) cease to exist beyond a certain finite time, where 
k(x) — > oo. This will be referred to as the finite time singularity. 

For given initial conditions we can determine the singularities as function 
of time and thus the asymptotics of u. We assume that u (x) > (the op- 
posite case will be treated in the next section), say monotonically increasing 
up to x = x* and monotonically decreasing for x > x*. An example is 

u (x) = xe~ x2 (36) 
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- a localized initial disturbance. Note that the initial conditions can always 
be scaled by a constant, say u = Av if time is also scaled as t = Tt such that 
TA = 1. 

The turning time t = a is found from 

k'(x) = (1 - 2f{x)t)-"/ 2 {\ + (3xf - 2f)t) = (37) 

In particular the time cii of the first turning point occurs at the minimal 
time for which 

1 + (3xf (x) - 2f(x))t = 1 + (3x 2 u' + xu )t = (38) 
and in the special case (|36| ) the minimum occurs when 2 — 8x 2 + 3x 4 = or 



Xl = V (4 + Vl0)/3 w 1.54 (39) 

(x = \J (A — \/IC))/3 must be discarded since t < 0) and 

1 

ai = ———7 — s r w 0.44 (40) 

2/(x 1 )(3x?-2) 1 ; 

As t increases we see from ( |58D that the turning point x has to occur in 
the regime where uq < and thus x > x*. For the case (|36|) we get as t — >■ oo 



we get Zoo = ^2/3 w 0.816. 

The finite time singularity occurs when /c(x) — > oo, which by (|35|) is 
equivalent to 

( = /5 = 57M (41) 

The smallest value (5\ occurs at the maximum of f(x). For the case (|36| ) this 
is x = 1 and fa = e/2 w 1.36. 

For very large times the solution makes sense up to 

x c (t) « (2t)- : / 2 (42) 

which is much smaller than Xoo > x*. At large times thus only the part 
< x < x c will contribute to the solution and the turning point becomes 
irrelevant. We can therefore use the form (|33|) expanded for small x to find 



u 



(k, t) = x c (tf 3 k-^ 3 w (2t)- 2 / 3 k-^ (43) 
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3.3.2 Shock conditions 



Due to the multiple valuedness this solution is cut off by a shock, where 
it jumps approximately to 0. The appropriate shock conditions have to be 
found by using the appropriate conservation law. It is reasonable to assume 
that this is the conservation of energy, which is the fundamental ingredient 
in the shell model. In terms of the continuous wavenumber, the conservation 
of energy is 

1/4^° < 44 > 

which can be written as the local conservation law 

= ~3k (45) 

i.e. 

v 2 

q-) t = -(2u 3 k) k (46) 

As long as no shock occur ( |i"6| ) and (^) are equivalent, but crossing a shock 
the appropriate shock condition is 

y = d i = j\ (47) 

dt [q] y ' 

where [.] denotes the discontinuity across the jump and V is the velocity of 
the shock. Using energy conservation gives 

" = = 2rf < 48 > 

3.3.3 Asymptotic forward pulse 

Combining the shock condition (fig) with the asymptotic form ( 43|) gives for 
the edge of the shock k e (t) 

^ = 2(2t)- 2 / 3 (A; e ) 5 / 3 (49) 

which leads to 

k e (t) = (A-ct 1 / 3 )- 3 / 2 (50) 

(where c = 4 x 2 _2//3 and A is an unknown constant). Thus the shock position 
diverges at the finite time t* = (A/c) 3 where the spectrum becomes of the 
Kolmogorov form all the way to the largest k and only decays as i~ 2 / 3 . 
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3.3.4 Backwards propagation 

In the backwards direction we use the fact that 

a = Reu = \/3lmu = \/3b (51) 

Then 

bt - Qk 2 bb k = 2kb 2 (52) 

and now the substitution u = —6/2 gives us back fl28| ) for u. Thus solving the 
propagation to smaller k can be done by taking a negative initial condition 
for fl28|) and we thus have the solution as fl3"4|) -fl35|) with u < 0, i.e. of the 
form 

uo(x) = —xe~ x (53) 

In this case the finite time singularity, where k(x) — > oo never occurs 
since 1 — 2ft > for all x, t > 0. The turning points, which create shocks 
are again found by k'(x) = 0. The analysis is very similar to the forward 
case, but now, for the initial condition (|53"D the first shock appears at x\ = 



(4 - V10)/3 w 0.53 with e*i w 2.04. 
For £ > «i the turning points occur at two distinct values of x, x a and 
Xb > x a . For t — > oo we see from ( P5D that they are solutions of 3x 2 -Uq + 
x«o = which means that x a approaches the left hand edge of the interval 
of support for u (which is zero for ([53])) and Xb < x* since u\xb) < (for 
(53|), — > ^2/3). (It can be noted in passing that the maximum for w also 



occurs very close to xj,, but at slightly larger x. The curve thus never crosses 
itself and only develops a cusp in the limit t — > oo.) 

3.3.5 Asymptotic left-going pulse 

Asymptotically the only relevant regimes of initial data are < x < x a and 
Xb < x < oo. For very small x we simply expand in x to get u(k). The 
interesting regime is Xf, < x < oo and here the main variation in u comes 
from the regime where —2ft ^> 1. There 

ukt = u (l - 2ft)-' i ' 2 x{l - 2ftffH = ft(l - 2ft)- 1 -> - X - (54) 



or 



(55) 
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Again, the shock condition ( f48|) gives (since u jumps to a very small value) 

§ = = (56) 

with the solution 

k e (t) = j (57) 



Asymptotically the solution thus looks like ( p5j) down to k = k e given by 
( |57D where it jumps to a small value (0(l/t)). Note that the value of u on 
the edge of the shock remains constant 

"« 53 ~h (58) 



4 Conclusion 

We have shown that an initial pulse on top of the zero solution of the inviscid 
Parisi equations splits into two part, one for each direction. The asymptotical 
solutions for the two pulses are: 

Uright = (2t)- 2 / 3 k- 1 ^ (59) 

uieft = {tk)~ down to k ~ C/t (60) 

both decaying in time. Thus a burst created on top of the zero solution not 
only travels down the inertial range until it is dissipated by viscosity (here, 
by the continuum), it also travels upwards, to the smallest k. The burst does 
not remain localized in fc-space, but is distributed over the /c-range, and then 
decays. It is interesting to note that the initial disturbance creates both a 
forward and an inverse cascade. 

We would like to thank Jakob Langgaard Nielsen for helpful discussions. 
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